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Any experiment attempting to verify the presence of entanglement in a physical system can only 
generate a finite amount of data. The statement that entanglement was present in the system can 
thus never be issued with certainty, requiring instead a statistical analysis of the data. Because 
entanglement plays a central role in the performance of quantum devices, it is crucial to make 
statistical claims in entanglement verification experiments that are reliable and have a clear inter- 
pretation. In this work, we apply recent results by M. Christandl and R. Renner [T] to construct a 
reliable entanglement verification procedure based on the concept of confidence regions. The state- 
ments made do not require the specification of a prior distribution, the assumption of independent 
measurements nor the assumption of an independent and identically distributed (i.i.d.) source of 
states. Moreover, we develop numerical tools that are necessary to employ this approach in practice, 
rendering the procedure ready to be applied to current experiments. We demonstrate this technique 
by analyzing the data of a photonic experiment generating two-photon states whose entanglement 
is verified with the use of an accessible nonlinear witness. 

PACS numbers: 03.67.-a, 03.65. Ud, 03.67.Mn 



I. INTRODUCTION 

Entanglement plays an essential role in various quan- 
tum information processing tasks [2HS] and experimen- 
tal verification of entanglement is crucial for testing and 
characterizing quantum devices such as sources and chan- 
nels As these devices move closer to the realm of 
practical technologies, our ability to perform reliable en- 
tanglement verification tests becomes increasingly impor- 
tant. Correspondingly, many theoretical and experimen- 
tal procedures for entanglement verification have been 
proposed (see [2l f5l and references therein) and their im- 
provement and development remains an active area of 
research [71[S]. 

Any entanglement verification procedure can be 
thought of as a series of measurements on a physical sys- 
tem followed by an analysis of the outcomes. The data 
obtained from these measurements is necessarily finite 
and therefore the claim that entanglement was present 
can never be issued with certainty. More precisely, there 
will always be a non-zero probability that the data was 
produced from a separable state, regardless of what the 
data may be. We are thus forced to provide statistical 
statements that quantify our confidence that entangle- 
ment was indeed present. Naturally, the procedure that 
leads to these statements should have a clear interpreta- 
tion, should not rely on unwarranted assumptions about 
state preparation and be readily implementable in prac- 
tice 0. 

The most widely used approach consists of computing 
the standard deviation of measured quantities and us- 
ing these as error bars to specify the uncertainty of the 
reported values [5j. However, there are several concep- 
tual issues with this approach [TOl [H] , including the fact 
that it can lead to counter-intuitive results (121 and is 



known to be inadequate to deal with nonlinear expres- 
sions [13]. This strongly asks for better alternatives and 
consequently other approaches have been recently sug- 
gested (see e.g. [11]). 

In this paper, we apply the work of Christandl and 
Renner on quantum state tomography [T] to formulate 
a reliable method for analyzing the data of entangle- 
ment verification experiments. As shown in Ref. [T], the 
method does not rely on the specification of a prior dis- 
tribution of prepared states nor on the assumption that 
they are independent and identically distributed. Addi- 
tionally, it is suitable for experiments performing arbi- 
trary quantum measurements and the final statements 
have a clear and well-defined operational interpretation. 
The approach relies on the concept of confidence regions: 
regions of state space that contain the true state with 
high probability [T]. 

Applying this method requires the specification of a 
region of state space for all possible measurement out- 
comes, an issue that is not dealt with directly in Ref. [T]. 
In our work we provide a recipe to assign confidence re- 
gions to data obtained from entanglement verification ex- 
periments that rely on entanglement witnesses. This as- 
signment requires the evaluation of a non-trivial inequal- 
ity for which we specifically develop numerical techniques 
to efficiently calculate it, rendering the entire method 
ready to be applied to current experiments. We demon- 
strate this fact by experimentally producing a family of 
entangled two-photon states whose entanglement is ver- 
ified by an accessible nonlinear witness (ANLW) [15) . 

The remainder of this paper is organized as follows. 
For the sake of completeness, we first briefly outline the 
framework introduced in Ref. [Ij and summarize some 
of its main results. We then proceed to illustrate the 
data analysis procedure that we build and elucidate the 
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numerical tools that we develop to perform the necessary 
calculations. Finally, we describe the experimental setup 
and analyze the results with our technique. 



II. CONFIDENCE REGIONS 

We now provide an overview of the main results of Ref. 
[1] and direct the reader to this work for further details. 
We begin by considering a collection of n + fc quantum 
systems described by a state p""'"'^, each system associ- 
ated with a Hilbert space H of dimension d. The mea- 
surement is performed only on the first n systems and is 
described by a general POVM consisting of a set {Bi} of 
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In the case of 



positive operators satisfying J2i Bi 
independent measurements of each of the systems, each 
element Bi will be a tensor product of n positive oper- 
ators acting on a single copy of the state. However, it 
must be clear that the formalism does not require this 
assumption: one should always think of this POVM as 
an arbitrary, generally collective measurement on H^". 
The role of the remaining k systems is purely operational: 
the goal of the entanglement verification procedure is to 
make predictions about the state of these remaining sys- 
tems. More precisely, we want to know if these systems 
belong to regions of state space that contain only entan- 
gled states. Note that n is the number of runs of the ex- 
periment, producing n systems which are then measured 
and the outcomes analyzed to build the predictions. 

Consider an experiment in which the predictions are 
made only for a subset of k' subsystems, k' < k. It was 
noted in Ref. J, that in the limit of A; — >■ oo, the reduced 
state of the n + k' subsystems — Tik-k' {p^^^) can 

always be described by a permutationally-invariant state 
of the form f P{a)a'^'-"'^''">da [H]. This corresponds to 
the usual independent and identically distributed (i.i.d) 
case in which many copies of a true state a are prepared 
according to some initial probability distribution P{<j). 
Thus, in the scenario of an experiment that can in prin- 
ciple be repeated an arbitrary number of times {k — ?> oo) 
and predictions are made on a sample of k' states, the 
above result in fact provides a justification of the i.i.d. 
assumption that is common in the literature. For conve- 
nience, we will adopt this point of view but remind the 
reader that the i.i.d. assumption is not necessary for the 
validity of the upcoming results ilj. 

The data analysis procedure we will employ is a map- 
ping that assigns a particular region of state space to ev- 
ery possible measurement outcome. Crucially, this map- 
ping must be specified before the experiment is carried 
out. The regions are deemed confidence regions if they 
contain the true state with high, user-specified probabil- 
ity. More precisely, for all i, denote by R{Bi) the region 
assigned to outcome Bi. This region will be a subset 
of the space of density matrices associated to V.. 

Then any prescribed region R{Bi) is deemed a confidence 
region with confidence level 1 — e if it satisfies the prop- 



Probs^ [a e R{B,)] > 1 - e, Vct, 



(1) 



where Probs- [a e R{Bi)] is the expected probability of 
success with respect to the distribution Tr (a^^Bi) of 
the measurement outcomes Bi. In this picture, statisti- 
cal statements take the following form: "We have applied 
a procedure that, with probability at least 1 — e, assigns a 
region containing the prepared state cr" . It is important 
to emphasize that this probability refers to the success of 
the procedure before any measurements are carried out: 
in the end, the original input state a is either definitely 
contained in the assigned region or not. The quantity 
1 — e should thus be interpreted as the confidence level 
of the statement that the state is contained in the as- 
signed region. This statement is valid for all possible 
states and outcomes and does not depend on extra as- 
sumptions about state preparation nor on the prior dis- 
tribution P{cr). This fact makes the procedure reliable 
and robust even in the cryptographic scenario in which 
a might have been chosen maliciously [1]. 

A main result of Ref. [T] was to provide a criteria 
to determine whether a given mapping from outcomes 
to regions succeeds in constructing confidence regions. 
This result is summarized as follows. Firstly, for each 
measurement outcome define the function 



where 
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Tr ((t'^"B, 
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(2) 



V(H) 



Ci{a)da 



is a normalization constant. The function Tr((T'*"i3i) is 
usually referred to as the likelihood function, so that /ij (cr) 
is simply its normalized version. Furthermore, let {F,} 
be a collection of subsets of 'D{'H), where the number of 
these regions is equal to the number of POVM elements 
{Bi}. For each set F^ define the enlarged set 

Ff = {cr : 3cr' e Ti such that F(cr,cr') > - ^2}^ (3) 
where F{a, a') = Tr ( \/~^foa\fo\ is the fidelity and 



In - + - l)lnn 



If for all possible outcomes Bi it holds that 



with 



IJ.i{a)da > 1 



(4) 



(5) 



(6) 



then the assigned regions F,^ are confidence regions with 
confidence level 1 — e (Corollary 1, \V). In equation ([5|, 
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da is the Hilbert-Schmidt measure: the flat measure on 
the set of density matrices of dimension d induced from 
the Haar measure on the set of pure states of dimension 
d X d [T7]. It must be noted that the polynomial factor 
2n^^ -i)/2 is an improvement on the term appearing 
inRef. p. 

The above condition ([5| can be more conveniently cast 
by referring directly to the quantity 1 — Jp ^i(a)da and 
making a direct comparison with the term e/cn.d- This 
can be achieved by instead integrating over the comple- 
ment regions Ti = {a : a ^ Ti}. Therefore we define 

e2{Bi,ri) := J_fi,,{c7)d(T 

_ /r-A(cr)c?o- 
Iv(n) A(cr)dcr 

For convenience, we will drop the explicit dependence 
on Bi and Ti from e2{Bi,Ti) whenever it is not neces- 
sary, while keeping in mind that its value will depend on 
the measurement outcome and the region assigned to it. 
Condition ^ can then be more conveniently cast as 

£2 • Cn,d < e- (8) 

In summary, the assigned regions {Fj} determine whether 
criteria Q is satisfied for some fixed value of e and when- 
ever it is, the enlarged regions F,^ constitute confidence 
regions. It is these latter regions that we assign to each 
individual outcome in our data analysis procedure. 

It is very important to note the role played by the 
polynomial factor Cn,d and the enlarging parameter S. 
Because the dimension of the Hilbert space d is fixed for 
a given experiment and typically large, the factor Cn^d 
will be a high-order polynomial in the number of runs 
n. Satisfying condition ([s]) will require £2 to be much 
smaller than the value of e that quantifies the confidence 
of the procedure. This can be problematic for small n 
but will play only a minor role for larger values because 
£2 decreases exponentially in n whenever the maximum 
of the function /i,; (cr) is contained in the region Fj [1. 

On the other hand, the size of the complement region 
Fj increases as S grows larger, implying that large values 
of 5 result in larger values of €2- In particular, whenever 
6 > 1 (which can occur for sufficiently low n) it will hold 
that the region F^ will be equal to the entire state space 
'D{H) and consequently £2 = 1. Thus, for a fixed confi- 
dence level, the value of n for which (5 = 1 sets a lower 
limit on the number of runs of the experiment that are 
required to verify the presence of entanglement. This is 
illustrated in Fig. [ij These features indicate that in this 
framework, it is usually necessary to accumulate large 
amounts of data in order to reliably report the presence 
of entanglement. 

We have in hand a method to verify whether a set of 
prescribed regions are in fact confidence regions. The 
question then remains of how to choose these regions in 
the first place, an issue that is not addressed in Ref. [T]. 
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FIG. 1: (Color online) Contour plot of 5 as a function of 
the confidence and number of runs n. The red region to the 
farmost left represents the case when S > 1, illustrating a 
lower bound on the number of runs that must be performed 
to achieve a certain value of e, quantified by the quantity 
— logjQ e. In practice, even larger values of n will be required 
to meet a desired confidence. 



Although the results of Christandl and Renner were orig- 
inally targeted at quantum state tomography, we will in- 
stead apply these results in the context of entanglement 
verification. We now describe a procedure for entangle- 
ment verification that fully specifies how to assign confi- 
dence regions in terms of entanglement witnesses. 



III. ENTANGLEMENT VERIFICATION 
PROCEDURE 

The goal of an entanglement verification experiment 
is to determine whether a prepared state is entangled or 
not with the highest possible certainty. In the language of 
confidence regions this translates to the task of deciding 
with the highest level of confidence possible whether the 
prepared state lies in a region consisting only of entan- 
gled states. Reconstructing the state of a general quan- 
tum system is experimentally demanding, as the number 
of required measurement settings will in general increase 
exponentially with the number of qubits Moreover, 
even if a given state is completely specified, deciding con- 
clusively whether it is entangled is computationally de- 
manding and it is in fact an NP-hard problem in terms 
of the dimension of the system [12]. A way to circum- 
vent these issues is to focus on entanglement witnesses, 
the use of which has become an increasingly popular tool 
both in theory and experiments [20Vl24j . thus playing a 
central role in the field of entanglement verification. 

A linear entanglement witness W is an observable sat- 
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FIG. 2: (Colour online) The region Vfy is fixed as the set 
of states detected by a linear entanglement witness W. This 
region can be seen as the set of states above the black line. 
Fixing implicitly defines a region Tw that determines 
if criteria ([8| is satisfied. This region is located above the 
dashed green line labelled Tw ■ The required numerical efforts 
are greatly simplified by realizing that the set of states F^ 
above the dashed red line constitute a subset of Fw as in 
Observation [T] 



isfying 

w{as) := Tr {agW) > for all cTj separable, 
w{ae) < for at least one entangled state cTe- 

Therefore, recording a negative expectation value is a 
conclusive indicator that the state must have been entan- 
gled. We refer to this as the state being detected by the 
witness. Calculating the expectation value of a witness 
operator can be performed efficiently for any state of ar- 
bitrary dimension. Moreover, experimentally determin- 
ing the expectation value of a witness generally requires 
considerably fewer measurement settings than a full re- 
construction of the state, making them very attractive in 
practical scenarios. 

One can also consider nonlinear entanglement wit- 
nesses [551 US] which can be viewed as powerful exten- 
sions of linear witnesses in the sense that they will always 
detect more states than their linear counterparts. Non- 
linear witnesses are described by their values w{a) which 
are nonlinear in the expectation value of the measured 
observables. They also satisfy the property that their 
value is negative only for entangled states. Moreover, 
accessible nonlinear witnesses were recently developed in 
[15j . demonstrating that their expectation value can be 
evaluated from the same data as the original linear wit- 
ness. Such nonlinear witnesses have also been recently 
applied in experiments |27) . 

The starting point of our procedure is the specifica- 
tion of an entanglement witness W and a POVM {B^} 
whose possible outcomes are sufficient to determine the 
expectation value of W. In our description w{a) refers 
to the value of a linear or nonlinear witness. Recall that 
in order to verify entanglement whenever it is present, 
we need to assign confidence regions that contain only 



entangled states. For this purpose, we define 

:= {a : wia) < 0} Vz (9) 

as the set of detected states. From the definition of an en- 
tanglement witness, F^ contains only entangled states. 
Our goal will be to report Tfy as the confidence region 
whenever possible. Going back to definition ([S]), notice 
that the set F^ is defined for a fixed F^. But in our pic- 
ture, we are interested in always reporting regions that 
contain only entangled states. Therefore, we alterna- 
tively choose to fix the reported region F^ and construct 
the smaller regions implicitly. From ([3| , it can be directly 
seen that if F^ is fixed, its corresponding subregion Fn/ 
is defined by 

Tw ■= W ■ max F{a,a') < v/l-<52}. (10) 

We are now ready to specify the mapping from outcomes 
to regions that constitutes the data analysis procedure 
for reliable entanglement verification. 

Data analysis procedure. To construct confidence 
regions with confidence level 1 — e in an entanglement 
verification experiment, apply the following rule to assign 
a region to each outcome Bf. 

1. Fix e. 

2. For each possible measurement outcome Bi, com- 
pute e2{Bi,rw) = JY^fJ.i{cr)da. 

3. If condition ^ holds, i.e. if €2 ■ c„^ < e, assign the 
set of detected states F^. Otherwise, assign the 
entire state space I?('H). 

Therefore, we assign only two possible regions: the set of 
detected states Tw or the entire state space 'D{'H). Note 
that the entire state space is trivially a confidence region 
for any given confidence level, so that our assignment 
indeed produces confidence regions. However, assigning 
the entire state space must be interpreted as the state- 
ment that for the given confidence level, it is not possible 
to certify that the set of detected states contains the true 
state. 

Even though the procedure is now completely specified, 
we are still faced with the difficulty of calculating €2- As 
a first step, we note that it is preferable to find a simpler 
way to characterize the set Tw- One way to do this is to 
find a subset of Tw that can be more easily described. 
We now show that such a subset can always be found 
in terms of a bound on the expectation value of a linear 
entanglement witness. 

Observation 1. Let W be an entanglement witness and 
let the number a > satisfy a > 2\\W\\ooS. Then the set 
Tq = : Tr(aW) < ~a} is a subset ofTyy- 

Proof: In order to prove the claim we only need to 
show that F'^{a,a') < 1 — S'^ whenever Tr (crW^) < —a 
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and Tr (cr'W) > 0. We begin by considering the following 
general inequality: 

\Tr[{a' -a)W]\ = ma' -a)\ 

< ||M/||oolk'-a||t,. 
<2\\W\\^^1-F^ia,a') (11) 

where we have used Holder's inequality 

|(^,W^)| < ||a||tr||l^||oo (12) 

and the Fuchs-van de Graaf inequality |28j 

\\a' -a\\t,< 2^1- F^a, a'). (13) 

Now let Ti{aW) = -a and Tr{a'W) = (3 for some 
a, (3 > 0. Inserting into (11) and rearranging we get 



J3+_a_ 



We want to find a condition on a such that F'^{a,(j') < 
1 — (5^ for any /3. This will occur whenever 



< 1 



13 + a 
a>2||H^||oo<5-/3. 



Since this inequality must hold for all /3, we can restrict 
ourselves to the worst case scenario of /3 = to obtain 



a>2||W^|U'^ 



(14) 



as desired. 



This result is illustrated in Fig. [2] Unfortunately, ob- 
taining a similar and useful result for nonlinear witnesses 
is difficult: the value of the nonlinear witness may differ 
greatly for two states even if their fidelity is high. 

Note that because Fq, C F^, it holds that 



/ Ci{(j)da > / Ci{a)d(j 



(15) 



since Ci{(T) > 0. Therefore if condition ^ is satisfied 
when integrating over F^, it will always be satisfied for 
the integral over Tw- 

Typically, it is possible to assign the set of detected 
states as a confidence region for very high confidence 
levels i.e. with e <C 1. Therefore, from now on we will 
quantify the confidence level of the procedure by the more 
appropriate quantity 



C 



(16) 



which we refer to as the confidence of the entanglement 
verification procedure. We further define this quantity to 
be zero whenever the assigned region is the entire state 
space ^{H). Thus, higher values of the confidence result 



in higher certainty that the state is contained in the set 
of detected states. 

From the description of the data analysis procedure, it 
should be clear that the crucial step is the computation of 
£2 : a highly non-trivial task that requires the normaliza- 
tion of the likelihood function as well as its integral over 
the implicitly defined set Tw In the following section, 
we construct and illustrate a series of tools developed to 
numerically evaluate an upper bound on 62, ensuring a 
method to verify condition (Isl). 



IV. NUMERICAL TOOLS 

There are several difficulties in calculating €2- An an- 
alytical approach is essentially intractable owing primar- 
ily to the high dimensionality of parameter space and the 
non-trivial geometry of the space of positive semi-definite 
operators [53]. Moreover, the region of integration Tw 
is not known in closed form but can only be cast as a 
black-box i.e. we can only ask whether a state lies in this 
region or not. Finally, we require any approximation of 
£2 to provide an upper bound on its value in order to 
ensure that the inequality £2 < c„.d£ is always satisfied. 

Fortunately, high-dimensional integration over black- 
box constraints can be handled with the use of Monte 
Carlo techniques. Most of these techniques are well sum- 
marized in |30j . In the Monte Carlo approach, the mean 
value of the integrand is approximated by the average 
value of samples randomly drawn from the region of in- 
tegration, which in conjunction with knowledge of the 
hyper- volume of the integration region can be used to cal- 
culate the value of the integral. Importantly, the number 
of samples can be chosen independently of the under- 
lying dimension and any constraint can be straightfor- 
wardly incorporated by checking whether a sample point 
lies within the constraint region. 

More specifically, the simplest version of a Monte Carlo 
technique to approximate a general integral of the form 
Ir fi^)^'^ involves a random sequence of N density op- 
erators {a I, CT2, . . . , (Tat} uniformly sampled inside R ac- 
cording to the measure da. By definition, the average 
(/) R of a function over a region R satisfies 



f{r7)d<J={f)R-VR 



(17) 



where Vr = da is the hyper- volume of the integra- 
tion region. The goal in Monte Carlo integration is to 
approximate the average of the function from the ran- 
dom sample. Namely, we approximate the value of the 
integral as 



fia)da 



1 ^ 



i=i 



Vr, 



(18) 



while keeping in mind that all sampled states lie in the 
integration region. Convergence to the true value of the 



6 



integral is guaranteed as iV — > cx) due to the law of 
large numbers |30j . A main drawback of this approach 
is that convergence can be extremely slow for highly- 
peaked functions such as >Ci(cr), since only very rarely 
will a state be drawn from the region surrounding the 
maximum of the function. This is particularly trouble- 
some for our purposes because an error in the calculation 
of 62 can lead to wrong conclusions about the confidence 
of the procedure. For this reason, we now introduce an 
approach that can be easily and efficiently implemented 
and provides an upper bound on £2- 

We first note that such a bound can be achieved by 
introducing a lower bound on the normalization constant 
J\f. Since the likelihood function is strictly positive, this 
can always be achieved by integrating over a subset R of 
V{'H), i.e. 



^ Jiy A(g)rf' 



la 



(19) 



We can use this fact to our advantage by restricting R 
to be a region around the maximum of Ci (tr) . Note that 
this maximum is unique and is in general achieved for a 
convex set of states [13]. Ideally, this region should be 
chosen to satisfy Jj^Ci{a)da « ^^^^^ Li{(j)da in order 
to provide a tight bound, but this is not necessary as 
the bound is guaranteed to hold for any R. Additionally, 
because the likelihood function is more flat around the 
maximum and R is much smaller than Vil-L), drawing 
random states within R will greatly improve the conver- 
gence of a Monte Carlo integration. 

We now illustrate how this region R can be constructed 
from a hyper-rectangle in parameter space. Following 
the convention of [31 , we begin by parametrizing any 
state (T e ^{T-L) in terms of the real- valued Bloch vector 
T = {ti,T2, . . . ,r<j2_i) as 



1 



(20) 



where the operators {Aj} are an orthogonal set of 
traceless Hermitian generators of SU{d) satisfying 

Tr ^Aj ^ =1. Any operator written in such a form 

is immediately Hermitian and of unit trace but may 
be non-positive for some vectors t. Thus, it will be 
important to keep in mind that not all possible vectors 
yield valid density matrices. With this parametrization 
the likelihood function will be a function of the Bloch 
vector Ci{a) = Ci{Ti,T2., ■ ■ ■ ,tj2_i). Our goal will be to 
define a region around the maximum that contains only 
valid states for which the value of the likelihood function 
is sufficiently large. 

Construction of integration regions. To construct 
a region R to be used in an approximation of the normal- 
ization of the likelihood function, perform the following: 



1. Calculate the maximum value of the likeli- 
hood function C™^^ and any vector r* = 



tained. 



. ,r^2_]^) for which this maximum is at- 



2. Find, for all j, the lowest possible quantities x~' > 
such that C,{tI,t^, . . . ,t* ± 2:=^, . . . , t*2_ J = 
£max 1^ for some fixed number 77 > 0. If no such 
values can be found for some j, let a;^ = 00. 

3. Find, for all j', the highest possible quantities > 

is still 



such that a{Ti , , . . 
a valid density matrix. 



4. Define — imii{x^ ,y^}. Then the integration 
region R is equal to all the valid density matrices 
within the hyper-rectangle r defined by r = {t : 



,Vj}. 



This construction is illustrated in Fig. [3j Note that 
the task of maximizing the likelihood function can be per- 
formed efficiently and is routine in the context of quan- 
tum state tomography. A good choice of rj will in general 
depend on each individual problem, but it should be cho- 
sen large enough to include only regions that contribute 
significantly to the integral. 

Once the hyper-rectangle has been constructed, it is 
straightforward to perform the Monte Carlo integration 
by sampling uniformly within the rectangle, while keep- 
ing only operators in that sample that are valid den- 
sity matrices. Let these sampled states form the set 
{fJi, (72: • ■ ■ J fTv}- The target integral is then given by 



Ci{a)da 



1 ^ 

Ci)R-VR. 



Vr 



(21) 



Because typical values of the likelihood function are ex- 
tremely small, it is preferable to work with the logarithm 
of the function and use the identity 

log {a + b)= log[exp(log a - log 6) + 1] + log b (22) 

to add the values of Ci{<jj) at each step of the algorithm 
and determine {Ci)R as in equation (21). 

In order to calculate Vr, we use the fact that the 
Hilbert-Schmidt metric on the space of quantum states 
generates the Hilbert-Schmidt measure [ST . The Hilbert- 
Schmidt distance between two density matrices is given 
by 

Dhs{(^1,(^2) = Ikl - 0-2||2 = VTr [((Ti - (72)^]- (23) 

This correspondence between metric and measure implies 
that the volume of the hyper-rectangle r can be found 
in the usual sense as the product of the length of its 
sides with respect to the Hilbert-Schmidt metric. More 
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specifically, let af = cr(Tj*, . 



,T*,_J. Then 



the length Ar^ of the jth side of r is given simply by 



= \Vt^j+rJ\,\\2=r+ +' 



(24) 



where we have used the fact that the operators \j are nor- 
malized with respect to the Hilbert-Schmidt inner prod- 
uct. The hyper- volume Vr of r is then given by 



n A., 



(25) 



This correspondence is also useful in generating a ran- 
dom sample, as one needs only to obtain a random num- 



ber within the intervals [t* 



Because not all 



operators in r are valid density matrices, Vr is in general 
larger than the hyper-volume Vr of the integration re- 
gion R. However, one can estimate R from knowledge of 
the fraction / of the randomly drawn operators that are 
valid density matrices. The relationship between these 
quantities is 



(26) 



which can finally be inserted in (21) to provide the nu- 



merical calculation of the target integral 



Ci{<j)d<j 



1 ^ 



/•JlArfc. (27) 




FIG. 3: (Color online) Construction of integration regions. 
We imagine a two-dimensional section of parameter space 
characterized by the variables ri and T2. Only the region 
inside the triangle contains valid density matrices and con- 
tours of Ci{(j) are shown in the background. To construct the 
integration region we do the following: 1. Find the maximum 
of the function and a state for which it occurs, in this case 
(Ti,r|). 2. From this maximum, find the displacements a;^ 
and such that the value of the function is decreased by a 
specified amount, in this case corresponding to the 6th con- 
tour line. 3. Find the displacements and that define 
the points where the boundary of valid states is met. 4. By 
choosing the minimum of these quantities in each direction, 
we construct a rectangle (dashed) and the integration region 
is the intersection of this rectangle with the space of valid 
density matrices. 



k=l 



To calculate /, it is sufficient to verify how many of 
the drawn operators are valid density operators and di- 
vide this number by the total number of randomly drawn 
operators. 

One could imagine that a similar technique could be 
used to calculate the integral J^r^ Ci{a)da appearing in 
the definition of £2- Unfortunately, this would greatly in- 
crease the computational efforts in the construction of R 
since one must additionally ensure that each of the drawn 
samples lie in Tw- Additionally, in this case restricting 
the integration region results in an incorrect lower bound 
on €2- Instead, we can construct an upper bound on this 
integral via the maximum of the likelihood function as 



< ( max A(a) ) (28) 

where Vxi{-h) is the Hilbert-Schmidt hyper-volume of the 
entire state space. This volume was calculated explicitly 
in |31| for Hilbert spaces of arbitrary dimension. We can 
then combine this result with our previous bound on the 



normalization constant to provide an overall upper bound 
on £2 ■ Since this value will be typically very small and in 
order to make a direct comparison with the confidence, 
we will henceforth refer to the logarithm of £2 for which 
we now have the inequality 



logio £2 < log 
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Vgiy A(g) Vv{-H} 
(A) 7? Vr 



(29) 



Of course, the average of the likelihood function over 
Tw will generally be much smaller than the maximum 
over this region, making the bound very loose. However, 
in practice this is not a problem because the above bound 
on £2 is dominated by the much larger differences between 
the global maximum of the function and its maximum 
over Tw More specifically, for experiments with a large 
number of runs (large n), it will typically hold that 



(30) 
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so that 



J 1 

niax^gr^£,(g: 




(31) 



and the value for log^o ^2 is not altered significantly by 
the loose bound. 

The final quantity we must be able to calculate is the 
maximum of the likelihood function over Tw- This again 
is a non-trivial global optimization problem involving a 
black-box constraint. As in the case of integration, the 
particular features of this problem impede the usual tech- 
niques and strongly ask for a Monte Carlo approach. To 
handle the optimization in the general case, we employ 
an adaptation to the quantum scenario of a simulated an- 
nealing algorithm (SA) based on the Metropolis-Hastings 
algorithm outlined in |32j . 

The SA algorithm is based on a biased random walk 
that preferentially moves to states with higher values of 
the objective function while still accepting moves to lower 
values with a probability governed by a global "tempera- 
ture" parameter. This last feature prevents the algorithm 
from being confined in local maxima. Unfortunately, this 
same feature makes the convergence slow, usually requir- 
ing many steps to reach close proximity to the maximum. 
For each step, one must additionally make the costly veri- 
fication that the states lie in the region of integration Tw, 
so it must be understood that run times are usually long. 
A detailed description of the algorithm is included in the 
Appendix. 

One drawback of the SA algorithm is that due to its 
stochastic nature, independent runs of the algorithm will 
generally yield different values. Moreover, by construc- 
tion these values cannot be larger than the global max- 
imum. In order to address this issue, one should esti- 
mate the numerical error by performing many indepen- 
dent runs of the algorithm and collecting statistics of the 
sample values. The usual choice is to calculate the stan- 
dard deviation of the values [30] and take this as the 
error. It is then important to ensure that condition ([s]) 
is satisfied well within this error. 

Nevertheless, we are still interested in obtaining a more 
efficient method to solve the maximization problem. We 
can achieve this for the case of linear witnesses by noting 
that for the subset Ta of Tw, it holds that 



max Ci{a) > max Ci{a) 



(32) 



since in that case Tw is a subset of Fq,. Therefore, we 
can provide a final expression for the bound on 62 as 



/ max^pj^ Ci{t 
logio £2 < logio TTTT 



[a) Vx>(-H) ' 



(33) 



where Fq, is defined as in Observation [T] This expres- 
sion has the enormous advantage that because the con- 
straint over Fq, is convex and Ci{a) is log-convex, the 
maximization of Ci{<7) over this region can be calculated 
with vastly greater efficiency using standard methods in 
convex optimization. 

We are additionally interested in reporting the highest 
possible confidence level, which corresponds to the case in 
which the equality €2 • — e holds. The value of e2 de- 
pends on the region Tw which in turn implicitly depends 
on e through the definition of the enlarging parameter (5, 
so that the above equality is in principle an equation to 
be solved for e. Unfortunately, there is no clear method 
of how to solve the equation directly, primarily because of 
the difficulty of calculating £2 itself. Instead, to achieve 
the highest possible confidence level, one must iteratively 
adapt the chosen value of e until £2 ■ c„^d ~ e while still 
satisfying the inequality (|8|. 

With these tools in hand it is now possible to apply the 
reliable entanglement verification procedure for both lin- 
ear and nonlinear witnesses. We now proceed to demon- 
strate the features of the method by applying the tech- 
nique to data obtained from an experiment generating a 
family of entangled two-photon states. The entanglement 
of these states is verified with the use of an ANLW. 



V. EXPERIMENT 

To apply our entanglement verification procedure 
to experimental data, we aimed to produce photon 
pairs in the maximally entangled states | $(</')) = 
^ {\HH) -h e*-^ IFF)), where \H) and \V) are defined re- 
spectively as polarization parallel and perpendicular to 
the optical table. A frequency doubled titanium-sapphire 
laser (80 MHz, 790 nm) was used to pump a pair of or- 
thogonally oriented 1 mm /3-Barium borate (BBO) crys- 
tals, as seen in Fig. |4] By pumping with diagonal po- 
larization \D) — {\H) + \ V)), the pump may produce 
photon pairs via type-I noncoUinear spontaneous para- 
metric down-conversion (SPDC) in either the first or sec- 
ond crystal |33j . Bismuth borate, a-BBO, and quartz 
crystals were used to ensure that each path was spatially 
and temporally indistinguishable, and the photon pairs 
were filtered using bandpass filters with a centre wave- 
length of 790 nm and a bandwidth FWHM of 3 nm. The 
single photon signal was measured with avalanche pho- 
todiodes (APDs) and coincidences were recorded within 
a 3 ns window. 

Single photons were detected at a rate of approxi- 
mately 200 kHz in each arm, with a coincidence rate of 
approximately 35 kHz when the measurements are set to 
HH or VV. A quarter-wave plate was tilted to introduce 
an arbitrary phase shift between horizontally and ver- 
tically polarized components, allowing control over the 
phase (jj- This setup constitutes part of the setup used 
for the experiment reported in j34j . The two-photon state 
was prepared for six values of 0, corresponding to a wave- 
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SPDC 




l^^l Polarizing beam splitte] 
^ Interference Filter 



Half-wave plate 
^ Quarter-wave plate 



FIG. 4: (Color online.) Experimental setup for producing 
!$(<;/))) = ^ {\HH) + e^'''\VV)) polarization states. Photon 
pairs are generated via type-I noncoUinear SPDC in a pair of 
orthogonally oriented BBO crystals and analyzed with wave 
plates and polarizing beamsplitters. The phase cf> is adjusted 
by tilting a quarter-wave plate. 



plate tilt range of twelve degrees and transforming the 
state from |$^) to |$+). 

Projective measurements were taken in three bases, 
corresponding to the eigenbases of the operators {ux 
fx, cr. 



XT ^ y 



Gz}- We will refer to the elements 
of these bases as \xi){xi\, \yi){yi\ and respec- 
tively. For example, the eigenbasis of Uz ® (Jz is given 

by = lifH), |Z2) = \iiv). \zz) - |Z4> = \yy). 

and similarly for the other bases. To verify the entangle- 
ment of these states, an accessible nonlinear witness was 
constructed from the linear witness VF = (l/4)(l-|-cra;(g) 
Ox — Oy®ay+(Jz®<Jz\ Following [T5j, the expectation 
value Woo (c) of the nonlinear witness for a state a can 
be expressed as 



(34) 



where 



c 

k : 



Tr[a(|V-)(V'-|f/)* 

Tr (cr[/*) 

Tr (aW) - ck, 



V2 



i\HV) 



\VH)) and the superscript t denotes 
partial transposition. By choosing U = (Jz ® cTz, this 
expectation value can be computed from the expectation 
value of the aforementioned operators and the nonlinear 
witness is accessible |15j . An accessible nonlinear witness 
was chosen because it detects these entangled states for 
most values of (j). 

In this experiment, all measurements are independent 
so that each element of the POVM {Bi} is a tensor prod- 
uct of the operators corresponding to possible individual 
outcomes. The likelihood function takes the form 



C,{a)=X{Tr{a\x,){x,\f'^ ■T,{a\y,){y,\r 
• Tr(cr|zj)(zj|)"'- , 




[ ■ ii\'' ■ vti ■ Vv ■ d'd ■ d'a ■ y\'d ■ a'a ■ il Lh n'L m 
a. (Jt. CTy 




1 



S PS I.IOtt 

(1) 




' H'H Vll Vv ■ D'D D'A A'D .VA LL Lh H'L Rli 

<j) ^ 1.177r 
(2) 



m 



[ IlV Vtl Vv D'D D'A A'D A'A LL Li; [i'L R'li 
0"z O",, 




(f) ^ 1.367r 

(3) 



(4) 





, ii n i jun i 

[ IlV Vtl Vv DD DA yVD AA LL RL RH 

0"Z <7,T <Jy 

cj) 1.727r 

(5) 



FIG. 5: (Color online.) Results of projective measurements 
on SIX states of the form ^ {\HH) + e''^|V"V)), corresponding 
to the eigenbases of the operators {a^ ® OxyO'y® (Jy^cr^® a^}. 
The approximate value of the phase is included for each case. 
Counts were integrated over 1 s per measurement setting. 



where is the number of times outcome |2:j)(a;j| is ob- 
tained and similar definitions hold for the other opera- 
tors, so that the total number of measurement outcomes 
is n = X]j=i "-^ + '^y + "■z ■ Note that in this case the mea- 
surement outcome is fully specified by the numbers 



{n^, n^, nl}. In the experiment, six states were prepared 
corresponding to six different values of the parameter cj). 
The measurement outcomes for each case are summarized 
in Fig. [5] 



(35) 



We have calculated the confidence as in equation ( 16 ) 
for the six preparations of the entire experiment. These 
results are illustrated in Table [Tj We can report very high 
confidences for almost all states, with the exception of 
state 4 for which condition ([8| cannot be satisfied for 
any value of e. This is not entirely surprising as this 
state presents the weakest correlations in the 
and bases leading to a value of the nonlinear 

witness that is closest to zero, as seen in Table [l] Thus, 
the outcomes for this case most closely resemble the ones 
that could be obtained from a separable state. This again 
is evidence that only large data which are clearly incon- 
sistent with separable states can lead to the reliable state- 
ments obtained from our procedure. 

Additionally, we are interested in understanding how 
the maximum achievable confidence depends on the total 
number of runs of an experiment. It is also important to 
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State 


Approximate phase 


Confidence 


Woo 




Total counts 


Confidence (Fq) 


Confidence (Tw) 


1 


l.lOvr 


5150 


-23.0 




1500 








2 


1.17n 


2050 


-15.2 




3000 


18 


24 


3 


l.SGn 


410 


-3.4 




6000 


165 


200 


4 


1.54n 





-0.3 




15000 


300 


315 


5 


1.72n 


1819 


-5.8 




30000 


660 


700 


6 


1.89n 


4980 


-13.6 




60000 


1378 


1500 



TABLE I: Calculation of the confidence and value of the non- 
linear witness for all prepared states in the experiment. The 
total number of counts obtained in each case was roughly 
35,000. 



TABLE II: Calculation of the confidence for samples of dif- 
ferent size from the outcomes of experiment 6 based on Fw 
and Fa. 



-10 



-15 



-'20 

-25 



T T N T 

/ \ 
/ v • 

/ \ 
/ \ 




20000 30000 40000 
Total counts 



FIG. 6: (Color online) Value of the nonlinear witness Woo{4>) 
for the six states prepared in the experiment (dots). The value 
of the nonlinear witness for the family of states = (1 — 
p)\^((j})){^{(l>)\ + jl withp = 1/42 is shown in the background 
(dashed). This curve is included only to illustrate the values 
of (j> for which it is difficult to verify entanglement and should 
not be interpreted as a fit to the data. The value of p was 
chosen to adjust the scaling to the recorded values. 



gain insight on the cost of using the bound of Observa- 
tion [T] for hnear witnesses. For this purpose, samples of 
different size were randomly selected from the outcomes 
of experiment (6) in Fig. |5] That is, from the entire set 
of observations in this experiment (shown in Fig. [s]), we 
randomly selected a subset of all the data and interpreted 
it as arising from an experiment with a fewer number of 
runs (counts). The confidence was calculated for each of 
them using both regions Tw and Fq,, this latter being 
possible because this state is also detected by the linear 
witness. The obtained values using these two different 
methods is portrayed in Fig. [7] and Table [TTj 

The results indicate that, as a percentage of the to- 
tal confidence, the loss introduced by considering Fq is 
small. It is also clear that a large number of runs are 
necessary in order to report a non-zero confidence, in 
accordance to our understanding of the role of the en- 
larging parameter S and the polynomial factor Cn.d as 
discussed in section HH To estimate the numerical error 
present in the SA algorithm, we performed 20 indepen- 
dent runs of the algorithm for the data of state 1 and 
found this numerical error to be 1.85%. In all calcula- 



FIG. 7: (Color online) Confidence for random samples of dif- 
ferent size, quantified by the total number of counts. The 
confidences were calculated for Fw (triangles) and F^ (dots). 
These results illustrate that the bound introduced by consid- 
ering the subset F^ is small and is not an impediment to reach 
a very large confidence. In the case of 1500 total counts, the 
confidence is zero, consistent with our understanding that a 
large number of outcomes are needed in order to reliably re- 
port entanglement with our technique. Moreover, the data 
shows that the confidence is roughly linear in the number of 
outcomes. 



tions it was ensured that condition ([8| was satisfied by at 
least ten times this numerical error. In the construction 
of the integration regions a value of 77 = 10^ was chosen 
for all cases. Finally, the CVX package for specifying and 
solving convex programs |35| was used to numerically cal- 
culate the global maximum of the likelihood function, as 
well as its maximum over F„. 



VI. CONCLUSION 

In this paper, we have applied the work of M. Chris- 
tandl and R. Renner in Ref. pP to the case of entangle- 
ment verification. Through the concept of confidence re- 
gions, we have provided a procedure to make reliable and 
efficient statistical statements quantifying the confidence 
level of having entanglement present in a physical system. 
These statements have a clear operational interpretation 
and do not require the specification of a prior distribu- 
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tion nor the assumption of independent measurements 
or i.i.d. sources. We have shown that this method can 
be appHed in practice by developing specific numerical 
tools designed to calculate all necessary quantities. For 
the particular case of experiments relying on linear en- 
tanglement witnesses, we have shown that the procedure 
can be implemented efhciently using only plain Monte 
Carlo integration and convex optimization methods. The 
procedure is ready to be applied to current experiments 
as we demonstrated by applying the technique to data 
obtained from an experiment generating entangled two- 
photon states. High confidence values can be achieved 
whenever the data is strongly inconsistent with a sepa- 
rable state and the number of measurement outcomes is 
large enough. Our results thus provide an illustration of 
the techniques that must be employed in current experi- 
ments in order to obtain clear and reliable claims. 

It is important to note that this work assumes that 
there are no systematic errors in the measurements per- 
formed. In any real experiment, there will always be 
discrepancy between the intended measurement and the 
one actually performed, no matter how small this discrep- 
ancy is. These systematic errors can in principle lead to 
incorrect statements and a method to incorporate it in 
the framework must be pursued. Numerical techniques 
also invariably involve errors and these should also be 
clearly incorporated in the framework. Future research 
may lead to improved algorithms. Finally, let us note 
that it is often desirable to quantify the amount of en- 
tanglement present as opposed to just verifying it. Our 
technique can in principle be applied to such cases by 
reporting regions that contain states with at least a cer- 
tain amount of entanglement. Future work can focus on 
including this feature into the procedure. 
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Appendix 

Here we fully describe the simulated annealing (SA) 
algorithm. The algorithm is based on a biased random 
walk in state space that preferentially selects states with 
a higher value of the likelihood function at each new step 
of the iteration. However, it also accepts jumps to states 




100 200 300 400 500 600 700 800 900 1000 



Step number in SA algorithm 



FIG. 8: (Color online) Three independent runs of the same 
simulated annealing algorithm for the data of experiment 1. 
Although all parameters are identical in each case, the output 
is slightly different in each case due to the stochastic nature 
of the algorithm. 

with lower values with a probability that depends on a 
global parameter T, usually referred to as the tempera- 
ture because of its similarity with the physical tempera- 
ture in the annealing process of metallurgy. 

Below is a full enumeration of all the steps of the 
algorithm to calculate the maximum value of the 
likelihood function £(cr) over the set T^. A graphical 
illustration of how the maximum value of the function 
is reached as the algorithm progresses is found in Fig. 
[8] The random walk here described is based upon 
the quantum adaptation of the Metropolis-Hastings 
algorithm depicted in [32]. 

Simulated annealing algorithm: 

1. Select an initial value Tq for the temperature T as 
well as for the "step size" A. 

2. Generate a d x d— dimensional random state 
according to the Haar measure, where d is the di- 
mension of the underlying Hilbert space H. Trace 
out one of the subsystems to obtain the state ctq. 
If (To S T\Y continue to the next step, repeat oth- 
erwise. 

3. Randomly choose a 2 x 2 Hermitian matrix Hki in 
the following way. Pick two integers fc, I randomly 
from the set {1, 2, . . . , d}. If fc < ^ ^ Hki = \k) {l\ + 
\l){k\, similarly if k > I ^ Hki = -i\k){l\ + i\l){k\ 
and finally if = / ^ Hki = \k){k\ - \k + l){k + 1| 
(set k + l = lifk = d). 

4. Pick a distance 6 by sampling from a Gaussian dis- 
tribution with mean and standard deviation A. 
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5. Compute the state = exp{iHkid)\tp) . Trace out 
one of the subsystems of \tp') to obtain the state Cg. 

6. If (7q ^ Tw, repeat steps 2 to 5, contmue otherwise. 

7. Evaluate the ratio i? = log (/:(ct^))/£(cto)). Ifi?>0 
(£(0-0) > £((To)), let (Ti — CTq. Otherwise, flip a coin 
with bias p = exp{-| log(/:(ao)) - log(/:(a(,))|/r}. 
If "1" is obtained (which happens with probability 
p), again let cri = ctq, otherwise (Ti = cto- 

8. Repeat steps 2-6 times to generate a set 
{ci, (72, • ■ • 7 cttv} corresponding to N steps of the 
random walk. For each step, adapt the tempera- 
ture via the cooling rule T(s) = Tq/s where s is 
the step of the walk. The maximum value of C{a) 
over this set is the output of the algorithm. 

The performance of the algorithm depends strongly on 
the value of A and this value must be adapted throughout 
each step of the walk in order to maintain a fixed average 
acceptance ratio, i.e. the fraction of times we jump to a 
new state. Various values for these ratios are suggested 
|36j . Similarly, the choice of initial temperature is crucial. 
Its role is to prevent the algorithm from being stuck in 
local maxima by allowing it to escape such cases in the 
initial stages of the algorithm. The temperature is then 



reduced to ensure that convergence to the maximum is 
attained. Therefore, the choice of initial temperature and 
cooling rule is essential and varies for different cases. In 
practice, they must be chosen for each particular problem 
based mostly on experience. 

Finally, in order to check whether a new state belongs 
in Tw , it is necessary to determine the maximum fidelity 
of this state with any state in this set. For this purpose, 
we exploit the fact that the fidelity function is concave 
in both its arguments and that the restriction p G 
is convex for both linear and nonlinear witnesses. These 
properties allow us to employ the highly efficient tools of 
convex optimization to solve the maximization problem. 
Concretely, for a given state cr, we verify membership in 
TvF by solving the problem 

maximize F(a, a') 
subject to cr' e F^ 

where a' must be forced to be a density operator. The 
state cr is a member of Tw if the solution to this problem 
is larger than \/l — S'^. In our case, the CVX package for 
specifying and solving convex programs [35] was used to 
numerically solve the problem. 
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